Multi-Bunch Solutions of Differential-Difference Equation for Traffic Flow 
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Newell- Whitham type car- following model with hyperbolic tangent optimal velocity function in a 
one-lane circuit has a finite set of the exact solutions for steady traveling wave, which expressed 
by elliptic theta function. Each solution of the set describes a density wave with definite number 
of car- bunches in the circuit. By the numerical simulation, we observe a transition process from 
a uniform fiow to the one-bunch analytic solution, which seems to be an attractor of the system. 
In the process, the system shows a series of cascade transitions visiting the configurations closely 
similar to the higher multi-bunch solutions in the set. 

PACS numbers: 02.30.Ks, 45.70.Vn, 47.54.+r, 05.45.Yv 



I. INTRODUCTION 

Traffic flow in a no-passing freeway has extensively 
study with the car-following models. Some of these mod- 
els can describe the spontaneous generation of density 
wave along the road, the traffic congestion, as the col- 
lective motion of cars. The optimal velocity (OV) model 
[0,^ reveals the features of the congested flow. In a cer- 
tain range of the mean headway, cars gradually compose 
bunch of cars, which moves backward with a constant 
velocity. The most of the studies are achieved by the nu- 
merical simulations under the cyclic boundary condition, 
which solve the equation of motion of the model; 



Xn{t) =a[V{AXn{t))-Xn{t)] 



(1) 



Xn and Axn denote the position of the n-th car and the 
headway of the car, Xn-i—x„, respectively. Drivers try to 
adjust their car's velocity to the optimal velocity V{Ax) 
according to the headway. The constant a parameterizes 
the sensitivity in the adjustment, which should be less 
than a critical value to generate the congestion. 

In the density wave of the congested flow, high density 
region, the bunch of cars, alternates with low density 
one. The cars run almost same headway and velocity in 
each of the regions, and in the interface of the regions 
the density changes rapidly exhibiting the characteristic 
kink-like shape. The shape of the kink on the interface 
seems to be determined only by the sensitivity a and the 
OV function V{Ax), irrespective of the initial condition 
or the mean headway. In the vicinity of the critical point, 
the interface can be approximately described by the kink 
solution of the modified Korteweg-de Vries equation . 
For some particular choices of the OV function, the exact 
solution which describes the interface has been obtained 
[Q,^. However, in the numerical simulations, the genera- 
tion process of the congestion displays much complicated 
aspects. A number of bunches with various length may 
arise in the circuit, and the fusion of bunches is often 
observed. One can predict neither the number of the 
bunches nor the lengths of them. 

On the other hand, the traditional car-following model 



which is described by the first order differential dif- 
ference equation. 



{t + T)^V{Axn{t)), 



(2) 



has been studied in the traffic engineering for a long time. 
The reaction time r represents the time lag which it takes 
the car to respond to the change of motion. Since the 
OV model (Q) can be considered as the truncated Tay- 
lor expansion, i„(i -|- t) ~ Xn{t) + TXn{t), of Eq. (||), 
these models rather resemble each other in the qualita- 
tive behavior, especially in the generation of the steady 
traveling wave. It is shown numerically that Eq. (|^) with 
a hyperbolic tangent OV function has the congested flow 
solutions quit similar to ones for the corresponding OV 
model 0]. 

Early in the 90 's, Whitham ||^ showed that the model 
with V{Ax) given by an exponential function has the 
pulse-like exact solution for steady traveling wave de- 
scribed by the elliptic functions. Although his choice 
of the function V{Ax) admits no existence of the car- 
bunching solution mentioned above, he pointed out a 
crucial relation between the time lag r and the propa- 
gating velocity of the traveling wave. (See Eq. (I).) We 
may call it Whitham condition. Recently j9j, it is demon- 
strated analytically that a class of the car-bunching solu- 
tions represented by the elliptic theta functions exactly 
satisfies Eq. (^) with the hyperbolic tangent OV func- 
tion. 



V{Ax) = (, + ??tanh 



Ax — p 
2a 



(3) 



The authors of Ref. report that the Whitham con- 
dition holds in the numerical simulation for much wider 
choice of the OV function. 

In this paper, we will investigate the structure of the 
class of the exact traveling wave solutions of Eq. (||). 
Under the cyclic boundary condition, each solutions cor- 
responds to a periodic density wave with definite length 
and number of bunches. Generally, some solutions differ- 
ent in the number of bunches are possible. They will be 
mentioned as the multi-bunch solutions. We also perform 



1 



the numerical simulation and observe the generation pro- 
cess of the density wave which shows a relaxation to one 
of the analytic solutions. 

In the next section, we investigate the parameter space 
of the exact solutions and find the finite set of the multi- 
bunch solutions. In section III, allowed parameters for 
the multi-bunch solution are determined. These param- 
eters will be used in the numerical simulation in section 
IV. The final section is devoted to summary and discus- 
sions. Some mathematical formula for the elliptic func- 
tions and the exact solutions are summarized in the ap- 
pendix. 



II. MULTI-BUNCH SOLUTIONS IN A CIRCUIT 



We begin by summarizing our previous results 
the exact solution with a width parameter S: 



Xn (t) ^Ct -nh +A In 



^0 i'^t-J-J-x+S^'i) 
^0 {vt - 



2A 



5,q 



of 



(4) 



where "doiv^ q) is the elliptic theta function and q, A, A, v, 
d, C and h are ansatz parameters which characterize the 
solution, q is modulus parameter of the theta function, h 
is the mean headway of N cars in the circuit with length 
L = Nh. The traffic flow expressed by (||) displays alter- 
nate appearance of high density region and low density 
region. Thus, the traffic in the circuit is divided into sev- 
eral number of "bunches" , which move backward. The 
cyclic boundary condition, Xn+N{i) = Xn{t) — Nh, and 
the periodicity of theta function, -doiv + 1, q) — 'do{v,q), 
imply that N/X must be an integer, which coincides with 
the number of bunches Ub ; 



A 



N 



( nb : integer ). 



(5) 



Then, the "wavelength" A is the number of cars within 
a successive pair of high and low density region, ap- 
proximately. The width parameter S, which ranges in 
< 2(5 < 1, determines the ratio of low density region 
in a wavelength. 2(5 A and (1 — 2^)A are the approximate 
number of cars in the low and high density regions, re- 
spectively. 

The exact solution satisfies the equation of motion 
with the OV function (^). As discussed in ||^, the 
traffic model with given time lag constant r can admit 
the exact solution (H) only if the Whitham condition 



2i/At = 1, 



(6) 



is met. Under the condition, we find out four relations 
between the ansatz parameters and the coefficients of the 
equation of motion, ^, rj, p and (7, as stated in the ap- 
pendix (B14)-(B17). Among the seven ansatz parame- 
ters, the mean headway h is already determined by the 
system size L and N . Furthermore, by the Whitham con- 
dition, we find that v is not an independent parameter 



but proportional to the inverse of A. The remaining five 
parameters (?, A, A, 5 and C are subject to the above 
four relations. It may seem that a parameter remains 
free and the system have a one-parameter family of the 
exact solutions. However, we notice that A is discretized 
as (^) so that the bunch number nb is an integer. More- 
over, the number of the bunches can not exceed the total 
number of the cars N at least. (We will see that much 
stronger restriction exists for the maximum number of 
the bunches later.) Consequently, the ansatz (^ gives us 
the finite set of the exact solutions which correspond to 
the nb-bunch states of traffic congestion. We call them 
the "multi-bunch solutions" . 

To examine the multi-bunch solutions, let us precisely 
analyze the four relations (B14)-(B17); 
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2t dp ""M^S-P)' 



Ap d 



2t dp 
p = h- A\ii 
a = A. 



In- 



d\{P) 



?9i(2^-h/3)i?i(2(5 
^i(2b-P) 
^i{26 + py 



(7) 

(8) 

(9) 
(10) 



Here and hereafter, A is changed to a much convenient 
variable, /?, which defined in the appendix as 



(11) 



Since the variable is also equal to nb/(2iV), it will be 
mentioned as the "bunch parameter". (Note that, in the 
expressions (|^) and (^) , /3/r is substituted for v by using 
the Whitham condition.) The equation ( |l(]|) simply says 
that A is identical to a. From the equation (|7|), we realize 
that the velocity parameter C can be expressed as an 
function of g, 5 and P (or A). Thus, the two nontrivial 
relations (||) and (||) should be solved to establish the 
existence of the multi-bunch solutions and to find out 
the allowed region of the parameters. 

First, we treat the equation (^), whose variables in 
the theta function can be decomposed by using the ad- 
dition formula as shown in the appendix. Replacing 
di M/d oiv) by the Jacobi's elliptic function \/k sn2Kv 
in ( B19 ), it becomes 



p_d_ 

'21(3 



In 
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sn^2KP sn^iKS 



(12) 



where A/rj (= a/rj) is denoted by Tc because of the fol- 
lowing reason. A linear analysis [p[pl|] gives the instabil- 
ity condition for a uniform fiow described by Xn'' (t) = 
V{h)t - nh, 



2TV'{h) EE r-sech^ 



2a 



> 



tt/N 
sin TT /N 



(13) 



(Note that the right hand side of (|T^) is almost equal 
to one, providing that N is not too small.) It follows 
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that Tj/a, the maximum value of 2V'(/i), determines the 
minimum value of time lag, Tc, for which the uniform 
flow becomes linearly unstable. One expects that when 
the condition, 



T > Tc = -, 

■n 



(14) 



is satisfied, there exist the uniform flow which decays to 
develop a congested flow. The Tc may be referred to as 
the critical time lag. By performing /3-derivative, the 
equation (|l^) can be easily solved with respect to sn^K5 
and gives 



Tc 2KI3 cn2Kl3 dn2A'/3 ' 



(15) 



1 - 



sn2K/3 



Since sii^AKS < 1, the right hand side of the above equa- 
tion should not exceed 1. Thus, we obtain the solvability 
condition, which guarantees the existence of S, as 



< 



sn2KP cn2Kl3 
2KI3An2KI3 ' 



(16) 



ranges in0</?<l/2, since nb < N. It can be easily 
checked that the allowed region actually disappears for 
T < Tc, which agrees with the linear analysis (14). The 
boundary curve in FIG. 1 crosses the /3-axis at /3o, the 
maximum value of (3. It is given by a solution of the 
equation obtained in the q ^ (fc ^ 0, A' ^ 7r/2) limit 
of the equality of (Eq): 



Tc _ -^.^ sn2A:/3o cn2A'/3o 
7 ^ 9™ 2KPq ATL2Kf3o 



sin 27r/3o 
2^/3o 



(17) 



Thus, /So is determined only by the value of Tc/t. The 
existence of the nontrivial upper bound of /3 implies that 
the number of the bunches is also restricted much more 
than rib < N . The maximum value of the number of the 
is given by 



(18) 



where [x] is the maximal integer that does not exceed 
X. Consequently, when the time lag r and the total car 
number N are given, the system has n^^^ exact multi- 
bunch solutions with the bunch parameters 
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FIG. 1. The allowed region of the solvability condition for 
Tc/r = 0.85869. The maximum value of bunch parameter 
Po is 0.149835. The vertical dashed lines indicate the ranges 
of q for some possible /3 with N = 20, which correspond to 
nb = 1, 2, 3, 4, 5. 
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2N' 



2N ' 



(19) 



When we specify one of the possible bunch numbers 
or /?'s, the solvability condition ( p^ ) gives us an allowed 
range of modulus parameter 



< g < g„ 



(20) 



where gmax is the value such that the equality of the solv- 
ability condition is held for the selected j3. The vertical 
dashed lines in FIG. 1 show the ranges of q for some 
possible (3 with N = 20, where n^^^ — 5. 



III. CONSTRUCTION OF THE MULTI-BUNCH 
SOLUTION 

In this section, we determine the width parameter 5 
and the modulus parameter q for any possible nb (or (3) 
to construct the rib-bunch solution. In the allowed range 
( pO| ) , the equation (^5|) has two branches of S as functions 
of q. One, which stays in < 2(5_ < 1/2, is given by 



2S- 



sn2Kf3 



Tc 2KP cn2KP dn2KP ' 



(21) 



1 - 



sn2K(3 



The inequality gives us the restriction for the number 
of bunches for given r. The right hand side of it only 
depends on /3 and q, where the complete elliptic integral 
of the first kind, K, and th e m odulus of Jacobi's elliptic 
functions, fc, are given in (A3) as the function of q. In 
FIG. 1, we show the allowed region of (|l|) in the (/3, q)- 
plane for Tc/t = 0.85869. Note here that (3{= nb/{2N)) 



where we take the branch < sn^^ < K for the inverse 
Jacobi's function sn~^. The other branch 1/2 < 2S+ < 1, 
which corresponds to K < sn^^ < 2K, is obtained as 
26+{q) = 1 — 2S^{q). These two branches are connected 
each other at q — gmax with 26 ^ 1/2 as depicted in FIG. 
2. On 25 — 1/2, the low and high density regions in the 
traffic flow contains same number of cars. Eq. (I) tells 
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us that such a "symmetric flow" can be yielded when the 
mean headway h is equal to p, the inflection point of the 
OV function, since t?i(l/2 + /3) = i?i(l/2 - (3) regardless 
of /3 or q. As long as 2S ^ 1/2, the traffic flow becomes 
asymmetric. For the first branch < 2(5„ < 1/2, the low 
density region contains more cars than high density one 
and the mean headway h is less than p. Contrary, for the 
second branch 1/2 < 25_|_ < 1, the congested region has 
more cars and h exceeds p. 




0.1 0.2 0.3 0.4 0.5 0.6 0.7 



FIG. 2. The width parameter 25±(q) for r^/r = 0.85869 
with A'^ = 20. Each curve corresponds to the possible bunch 
number rib in FIG. 1. 

By plugging 2(5±((7) into (||), we obtain the relation h 
and between the modulus parameter g; 



h±{.(]) = p + cr In 



79i(2(5±(g)~/3,g) 
^9i(2(5±(q)+/3,g)' 



(22) 



where /i_(< p) and /i+(> p) correspond to the two 
branches of 8. (Note here that /i_ (g) + h+{q) = 2p.) The 
combined entire function h{q) is shown in FIG. 3. We 
can determine the modulus parameter q, within a certain 
range of mean headway L/N, by solving h{q) — L/N, 
and obtain a exact multi-bunch solution of ^) by cal- 
culating S and other parameters through the equations 
above. 

In some case, we may find several solutions for q, which 
correspond to same or different bunch number. Although 
the exact solution can construct for each q, the stability 
of these solutions is other problem. We will discuss the 



issue in another work []12 
dicate the linearly unsta 



headway of the uniform flow Xn (<); 
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(0), 



|/i — p| < 2(TArccosh< 



/ T sin7r/A^ 



(23) 



Within the range, we may expect that the unstable uni- 
form flow will grow into one of the analytic solutions. 
However, the generation process of the multi-bunch so- 
lution from a given initial configuration and how one of 



the possible analytic solutions is selected may be more 
complicated problem. 
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FIG. 3. The mean headway h±{q) for Tc/t = 0.85869, 
p — 2, a = 1/2 with A'' = 20. Each curve corresponds to 
the possible bunch number Uh in FIG. 1 and 2. The horizon- 
tal dashed lines indicate the critical values of headway for the 
linear instability, which coincide with the q limits of h{q) 
for rih — 1. 



IV. NUMERICAL SIMULATIONS 

To investigate the generation process of the multi- 
bunch solution from a uniform flow, we perform numeri- 
cal simulations by solving the differential-difference equa- 
tion (S). We adopt an OV function 



V{Ax) = tanh(Aa; - 2) + tanh2. 



(24) 



where the coefficients are chosen as ^ = tanh2, 77 = 1, 
p = 2 and cr = 0.5 (r^ = 0.5). The time lag r = 0.58828 
gives Tc/t = 0.85869 < 1. We prepare 20 cars and ar- 
range them so that they form a uniform flow with a com- 
mon headway h = 1.88571 and an initial velocity V{h) = 
0.850233 in a circuit whose length is L = 37.7142. The 
uniform flow is sustained for a duration of t to prepare 
the initial function of the differential-difference equation. 
The solvability condition (|l^) tells us that maximally 
five-bunch mode is allowed. The modulus parameters 
q for the possible bunch numbers rib = 1, 2, 3, 4, 5 
are 0.70792140328755, 0.50113376, 0.3536167, 0.2418044, 
0.140292, respectively. Since the uniform flow is linearly 
unstable, it starts developing the density wave. In this 
case, the one-bunch mode is generated after enough re- 
laxation time, i ~ 6 X 10"*, ultimately. The result, which 
is displayed in FIG. 4 by dots, agrees quite well with the 
analytic one-bunch solution with q = 0.70792140328755 
as shown by the thin line in the figure. 



4 



Ax{t) 



2 . 8 




1.2 



10 20 30 40 

FIG. 4. The result of a numerical simulation with A'' = 20, 
Tc/t = 0.85869, V{A.x) = tanh(Aa; - 2) + tanh2 after Dots 
show the headway of a car. The line shows the corresponding 
one-bunch solution with q = 0.70792140328755. 

We observe much interesting phenomena in the gen- 
eration process of the analytic exact solution. At first, 
until t ~ 300, the initial uniform flow grows into a three- 
bunch configuration which closely resembles the exact 
three-bunch solution. Although the configuration keeps 
its shape for a while, it is gradually distorting and one 
of the bunches get closer to the other. At t ~ 4680, the 
fusion of bunches occurs and it is transformed into a two- 
bunch configuration. The two-bunch configuration lives 
about ten times longer than the three-bunch one. Sub- 
sequently, a bunch starts shrinking. The bunch is ab- 
sorbed by the other at t ~ 51560 as shown in Fig. 5, and 
then the exact one-bunch solution is finally accomplished. 
The simulation has been continued until t = 2 x 10^, 
at which the one-bunch configuration has still survived. 
These phenomena suggest that the one-bunch solution is 
an attractor of the system jl^ . 
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FIG. 5. The fusion process at t ~ 51560. A bunch starting 
to shrink at t ~ 48000 is absorbed by the other 



V. SUMMARY AND DISCUSSIONS 

We find a finite set of exact multi-bunch solutions for 
a car-following model with driver's reaction time lag in 
a circuit, which is described by a differential-difference 
equation under the cyclic boundary condition. When the 
circuit length and the total number of cars are given, we 
can calculate the possible number of car-bunches in the 
circuit and the profile of the density wave which describes 
the car-bunching. The numerical simulation shows that 
the linearly unstable uniform flow develops into the one- 
bunch solution, which seems to be an attractor of the 
system. In the relaxation process from the uniform flow 
to the one-bunch solution, the system stays for a while on 
the conflgurations which are similar to the higher multi- 
bunch solutions. For smaller number of bunches, the du- 
rations of these quasi stable configurations become longer 
exponentially (See Fig. 1 in p3[|). 

The series of cascade transitions among the multi- 
bunch configurations suggests that the multi-bunch so- 
lutions correspond to heteroclinic points of the system. 
However, any flow out of the one-bunch solution does 
not observed up to the present. As the other possibil- 
ity, each of the multi-bunch solution may be a kind of 
the Milnor attractor , which is unstable for any small 
perturbations, but globally attracts the orbits. To see the 
problem, we are going to have to investigate the stability 
and the attracting domain of the multi-bunch solutions 
more precisely. 

The density wave formation from a uniform flow on 
the OV model is very similar to the present model. The 
car-bunching and the fusion of bunches are observed also 
in the OV model. It may be possible to expect that the 
existence of the multi-bunch solutions and the qualitative 
feature of the cascade transitions are common to both 
models, although the analytic solutions for the OV model 
have not been obtained. 
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APPENDIX A: DEFINITIONS AND FORMULA 
OF ELLIPTIC FUNCTIONS 

The definitions of elliptic functions and some mathe- 
matical formula used in this paper are listed. The elliptic 
theta functions are defined as the infinite products; 
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Mv, 9) = 90 n (1 - 29'""' cos 2ttv + 

00 

q) = 2(7^/'' go sinTrw J]^ (1 - 2^^" cos27rw + q"^"), 

n=\ 

oc 

where <Zo = H " (^1) 

n=l 

The addition formula of the theta functions are 

+ w)dQ{v - w)'dl{Q) = ^l{v)dl(w) - d\(v)'d\{w), 

(A2) 

Jacobi's ehiptic functions and their modulus fc and com- 
plementary modulus k' can be expressed by the theta 
functions: 



sn2ii:w = ^ ^5^44, cn2ii:w \/l - sv?2Kv, (A3) 



dn2i^w = \J\ - kHv?2Kv, 
where K is the complete elliptic integral of the first kind. 

APPENDIX B: THETA FUNCTION FORMALISM 

Let us solve the equation of motion, 

Ax„(i) - p 



init + t) ~ S, + 1] tanh 
by the functional ansatz, 

Xn{t) =Ct-nh + A In 
using the variable notations, 



2a 



^^2A- 



(Bl) 



(B2) 



(B3) 



By virtue of the Whitham condition vt — /?, the velocity 
with time lag r, + r), can be written as 



Xn[t + T) = C + Av— In — — , 

dv Uq[v ~ 0) 



(B4) 



where we replace the time derivative d/dt with ud/dv. 
Converting the -y-derivative into d/ d5 again, we can apply 
the addition formula (A2) to the expression; 



d_ + _ ± 

dv "doi^ — dS 



ln[Mv + S)Mv - S)] 



-In [dl{v)4{S)-^l{v)dliS)] 



It allows the velocity to be a rational expression of '0q{v) 
and 

init + r) + Av 



We find that the headway, 



Axn{t) ^ h + A In 



Mv + P + S) Mv-(3-5) 



(B6) 



(B7) 



can also be expressed as a rational form of ^q{v) and 



„2X„ 



where 



Axn{t) — h 



2A 



(B8) 



(B9) 



By eliminating and 'd\{v) from (|B6|) and (|B8[), it 

will be shown that the velocity Xn{t + T) is equal to a first 
order rational expression of e^^" , which can be rewritten 
as an hyperbolic tangent function of the headway Ax„ {t). 

Performing the elimination, we get a differential- 
difference equation which the ansatz ([B2D satisfies; 



±n{t + T) = C + Av 



2X„ 



(BIO) 



where 

7v± = ± [^lis ± f3)idi{5)y - ^lis ± (3M{s)y] . (bii) 

We can transform N± into the expressions in terms of 
/3-derivatives, instead of S ones; 



±^^(0)^i(/3)i?i(2<5±/3)— In 



MP) ' 



(B12) 



(B5) 



As mentioned above, the differential-difference equation 
can be written by using the hyperbolic tangent function 
as 

.„(t + r) = C+-^- + ^ 

Aiy f N+\ / ID\^ 

+ -(^-^)tanh(x„ + -ln-j. (B13) 

Thus, we find the equations to which the ansatz param- 
eters are subject, by comparing the above equation with 
the equation of motion (Bl): 
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Av a Vi{lo + p) 


(B14) 






V = 


Aiy 
2 


(i/3 ??i(2d + /3)??i(2d - /?) 


(B15) 


p = 


h- 




(B16) 


a — 


A. 




(B17) 



In the first two equations, we can decompose the variables 
b, (3 in the theta functions by the addition formula; 



Av d 
~d6 



Av d , 
+ — 37 In 

4 d5 

Av d , 

7? = in 

' 2 dp 



dl{25) 



dl{25) 



(B18) 
(B19) 



dm di{25)_ 

For the first one, d/d[3 is converted into d/d5, to apply 
the addition formula, before the decomposition. All of 
the equations (B14)-(B17) can be shown that they are 
equivalent to the ones which we found in the previous 
paper although they are quite different in the ap- 
pearances. 
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